Viscosity calculated in simulations of strongly-coupled dusty plasmas with gas friction 
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A two-dimensional strongly-coupled dusty plasma is modeled using Langevin and frictionless 
molecular dynamical simulations. The static viscosity rj and the wave-number-dependent viscosity 
ri{k) are calculated from the microscopic shear in the random motion of particles. A recently de- 
veloped method of calculating the wave-number-dependent viscosity ri{k) is validated by comparing 
the results of ri(k) from the two simulations. It is also verified that the Green-Kubo relation can 
still yield an accurate measure of the static viscosity r; in the presence of a modest level of friction 
as in dusty plasma experiments. 



PACS numbers: 52.27.Lw, 52.27.Gr, 66.20.-d, 83.60.Bc 



I. INTRODUCTION 

Strongly-coupled plasma is a collection of free charged 
particles where the Coulomb interaction with nearest 
neighbors is so strong that particles do not easily move 
past one another. A widely used criterion to determine 
whether a plasma is strongly coupled is F > 1 [E 
where F is defined as the ratio of the potential energy be- 
tween neighboring particles and the kinetic energy. When 
F > 1, particles move slowly and are trapped by a cage 
consisting of a few nearby particles. If they escape the 
cages gradually, particles in a strongly-coupled plasma 
can flow, much like a liquid [s!]. However, if F ^ 10^, 
nearby particles that form a cage move so little that a 
particle inside the cage can seldom escape the cage; this 
condition is like molecules in a solid [1, Q • If a shearing 
stress is applied, cages in a solid are elastically deformed 
but can restore to their previous state, whereas cages in 
a liquid are disrupted and a viscous flow can develop. 

One type of strongly-coupled plasma is dusty plasma 
formed in the laboratory. A dusty plasma consists of four 
constituents: micron-size particles of solid matter (dust 
particles), electrons, ions, and neutral gas atoms [6|-l8|. 
The dust particles are strongly coupled amongst them- 
selves due to a large interparticle potential energy pro- 
vided by a large particle charge [ol. llOf. Several schemes 
have been used to confine charged dust particles us- 
ing natural electric field inside a plasma. One of these 
schemes makes use of a radio- frequency plasma [ll|, , 
with a horizontal electrode that provides a sheath elec- 
tric field that can confine and levitate dust particles in 
a cloud with only a few horizontal layers. If experi- 
menters introduce only a limited number of dust par- 
ticles, they can settle into just a single layer In these 
single-layer clouds, dust particles have negligible verti- 
cal motion, so that the cloud of dust particles is often 
described as a two-dimensional (2D) system [ol. Iisl- 16|. 
In this 2D cloud, the interaction between dust particles 
is a repulsive Yukawa potential Due to the large 

length scale and the slow time scales dusty plas- 
mas allow video microscopy to track individual particle 



motion flq. In dusty plasma experiments, elasticity in 
solids [14] and viscosity in liquids 16] has been observed 
and studied. However, strongly-coupled plasmas cannot 
always be classified as purely elastic or purely viscous. 

Dust particles experience several forces in the ex- 
periments. The electric force provides strong coupling 
amongst the dust particles as well as the levitation and 
confinement. Gas friction, due to dust particles moving 
relative to the rarefied gas, is the primary energy loss 
mechanism. The gas is usually so rarefied that it rep- 
resents only a small portion of the mass of the dusty 
plasma. Gas represents < 10% of the mass of dust 
in a 3D dusty plasma experiment at 400 mTorr [l9|, 
while 2D experiments have even less gas, with a pres- 
sure < 20 mTorr SHIIH. There is an ion drag force 
due to a steady flow of ions, arising from the same dc 
electric fields that provide levitation and confinement of 
dust particles. This ion drag force is parallel to the ion 
flow. Finally, in some experiments, laser radiation pres- 
sure forces are used to accelerate dust particles, for exam- 
ple to create macroscopic flows [13, 3, 21, 22] or simply 
raise the kinetic temperature of the dust particles with- 
out causing a macroscopic flow [3, 0, 3, 24 ] . This kind 
of laser heating method is one of several ways that exper- 
imenters can control F so that the cloud of dust particles 
behaves like a liquid or solid 13|, 1241427 ) . 

We assume that the Coulomb interaction amongst 
charged dust particles is the dominant mechanism for vis- 
cosity in laboratory dusty plasma experiments. Viscous 
transport of momentum occurs when the dust particles 
moving relative to one another in a shearing motion col- 
lide, causing some of their momentum to be transferred 
across the flow. We expect that collisions involving gas 
atoms will contribute less to the viscosity. Although the 
force of gas friction is effective in diminishing momen- 
tum of dust particles in the direction of their motion, 
there are two reasons it has little effect in transferring 
momentum across a flow of dust particles. First, the 
gas is rarefied so that it can carry much less momentum 
than a viscous solvent in a colloidal suspension (isl [29| , 
for example. Second, in a 2D experiment like l20|, a gas 
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atom that is struck by a dust particle is usually knocked 
into a direction out of the dust layer, so that there is lit- 
tle opportunity for a dust particle to push another dust 
particle indirectly through collisions with a gas atom. 

Here, we will refer to the viscosity as the static viscosity 
rj to distinguish it from viscoelasticity. In the literature of 
dusty plasmas, the static viscosity rj has been measured 
experimentally 16 1 and quantified in simulations 30l-33[. 
There are two ways to quantify the static viscosity. If 
there is a macroscopic velocity shear, the static viscosity 
can be calculated from the velocity flow profile 12, ll6, 
32 . [3^. On the other hand, if there is no macroscopic 
velocity shear, the microscopic shear associated with the 
random motion of particles can be used to calculate the 
static viscosity using the Green-Kubo relation 3l|, |33, 
111- 

Viscoelasticity is a property of materials that exhibit 
bot h liq uid-like viscous and solid-like elastic characteris- 
tics [35|. Most materials in reality are viscoelastic, such 
as wood, synthetic polymers, and human tissue [35|- Vis- 
cous effects correspond to energy dissipation, while elas- 
tic effects corresponds to energy storage. In general, liq- 
uids exhibit mostly viscous effects at large spatial and 
temporal scales, but they exhibit some elastic effects at 
small spatial and temporal scales [36l |. 

To characterize viscoelasticity quantitatively, it is com- 
mon to use either the frequency-dependent viscosity r]{uj) 
or the wave- number-dependent viscosity r](k). The latter 
characterizes materials at different length scales, and was 
introduced by theorists performing simulations 37-4ol|. 
The static viscosity rj is the hydrodynamic limit of the 
wave-number-dependent viscosity ri(k) when fc — > 0. In 
considering this limit, the relevant characteristic length 
scale for k is the interparticle distance, which is often 
measured as the lattice constant 6 of a perfect crystal. 

Viscoelasticity of stro ngly -coupled plasmas has been 
studied theoretically 41-43| and experimentally [20I. H. 
I45I I. The few experiments that have been reported for 
viscoelasticity of dusty plasma include a descriptive pre- 
sentation [44| and a characterization using a correlation 
function of the microscopic motion of dust particles (45| . 

In our recent 2D experiment (20| . a single horizontal 
layer of electrically charged dust particles was levitated 
in a glow-discharge plasma. The kinetic tem per ature of 
the dust cloud was raised by laser heating [23,123. View- 
ing from above, we recorded movies of particle motion, 
then calculated particle positions and tracked them to 
calculate their velocities. Based on the trajectories of 
particles, the wave-number-dependent viscosity r]{k) of 
the 2D dusty plasma was quantified using an expression 
we derived that accounts for gas friction. 

In simulations, the viscoelasticity of both 2D [20'| and 
3D strongly-coupled plasmas [13] have been studied re- 
cently. In this paper, we carry out further simulations 
for two purposes: to validate the ri{k) calculation method 
taking into account gas friction, as presented in (20j , and 



to assess the accuracy of the Green-Kubo relation for 
dusty plasmas with a modest level of gas friction. 

Simulations of strongly-coupled plasmas usually use 
the molecular dynamical (MD) method [13, 113] ■ Each 
particle is tracked individually, unlike the case of particle- 
in-cell (PIC) simulations, where aggregations of particles 
are simulated by a hypothetical super-particle [7,] . Track- 
ing individual particles is suitable because otherwise the 
dominant effects of strong particle-particle Coulomb in- 
teractions would be lost. Another difference is that in 
MD simulations, as compared to PIC simulations, Pois- 
son's equation is not solved. The only equation that is 
solved is the equation of motion for each particle, which 
is integrated to track particle trajectories. The result of 
the MD simulation is a record of all particle positions 
and velocities, which is the same kind of data that are 
produced in dusty plasma experiments. The interpar- 
ticle interaction that is assumed in MD simulations of 
strongly-coupled dusty plasmas is a repulsive Yukawa po- 
tential E3, 



<Pi.j = (3^(47reo?'ij) ^exp(-r,;j/AD), 



(1) 



where Q is the charge on dust particles. A/? is the Debye 
length, and Vij is the distance between the ith and jth 
particles. 

We list here additional parameters for the dusty 
plasma cloud. Because the dust cloud is 2D, we use 
an areal number density n and an areal mass density 
p = mn for the cloud, where m is the dust particle 
mass. We note that while the units for mass density 
and viscosity are different in 2D and 3D, the units are 
the same for the kinematic viscosity [l^, rj/p. Distances 
between dust particles are characterized by both the lat- 
tice constant b for a crystal or the 2D Wigner-Seitz ra- 
dius a = (n7r)~^/^ [i^. Time scales for collective mo- 
tion are characterized by the nominal 2D dusty plasma 
frequency ujpd = (Q^/27reoma'^)^/^. Gas friction is 
characterized by the damping rate Vf, which is the ratio 
of the gas friction force and the dust particle's momen- 
tum. 

We will discuss how to calculate 77 and ri{k) in Sec. II. 
In Sec. Ill, we will discuss our two MD simulation meth- 
ods, Langevin and frictionless. In Sec. IV, we will re- 
port new simulation data for ri{k) of 2D strongly-coupled 
dusty plasmas. We will validate our analysis method [2^ 
for calculating ri[k) in 2D strongly-coupled plasma with 
gas friction. We will also test the accuracy of the Green- 
Kubo relation with a modest level of gas friction as in 
our experiment. 
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II. METHODS FOR CHARACTERIZING 
VISCOSITY 

A. Static viscosity rj 

The Green-Kubo relation is widely used for calculating 
the static viscosity rj, based on the random motion of par- 
ticles. This method is used when there is no macroscopic 
velocity shear. The Green-Kubo approach assumes linear 
microscopic fluctuations and equilibrium fields in the sys- 
tem (sij . The assumptions of this approach are similar 
to those for the fluctuation-dissipation theorem 47, 48| . 
Previously, the Green-Kubo relation was generally used 
with data from frictionless simulations 30l 31. 33l sl. 43|. 
To calculate the static viscosity, first we calculate the 
stress autocorrelation function (SACF) 



Viy{t) for the ith particle. These are used to calculate the 



where Pxy{t) is the shearing stress 
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(2) 



(3) 



where i and j are indices for different particles, N is the 
total number of particles of mass m, = {xi,yi) is the 



position of particle i, Xi 



•^i Xj ; l/ij — Vi Vj : '^ij 



r.; — Tjl, and 



) is the interparticle potential. Second, 
we calculate the static viscosity rj from the Green-Kubo 
relation [31|], 



1 



VkT 



C^{t)dt. 



(4) 



Here, V is the simulation volume, which is replaced by 
the area of the simulation box for 2D simulations like 
those reported here. 

The Green-Kubo relation, Eq. (jH), is intended for use 
in equilibrium systems, but in this paper we will assess 
whether it can also be used in systems with a modest 
level of gas friction as in our experiment [l^. The dust 
particles in an experiment experience gas friction, in ad- 
dition to collisions amongst themselves, whereas only the 
latter are modeled in the Green-Kubo relation. We will 
carry out simulations, with and without friction, and ver- 
ify that Eq. (|4]) yields the same result in both cases. 



B. Wave-number-dependent viscosity rjik) 

The wave-number-dependent viscosity rj{k) character- 
izes viscous effects at different length scales. A method of 
calculating rj{k) from the trajectories of random motion 
of molecules in liquids has been developed In cal- 

culating rj{k) using this method, one starts with particle 
trajectories, such as Xi{t) and the perpendicular velocity 



transverse current, jy(k,t) = 'Y^f^^Viy{t) e'XY>[ikxi{t)]. 
The normalized transverse current autocorrelation func- 
tion [ii, EOl (TCAF) is then calculated as 

CT{k,t) = o;(fc,o)j,(fc,t))/(j;(fc,o)j,(fc,o)), (5) 

where the wave vector k is parallel to the x axis. (Here, 
k serves only as a Fourier transform variable, and is 
not intended to characterize any waves.) The wave- 
number-dependent viscosity of frictionless systems is cal- 
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usmg 



r]{k)lp = l/($fc2 



(6) 



where $ is a time integral representing the area under 
the TCAF after normalizing the TCAF to have a value of 
unity at < = 0. Generally, 77 (fc) diminishes gradually as k 
increases, meaning that viscous effects gradually diminish 
at shorter length scales. 

In [i^l, we generalized this expression as 



Ti{k)/p = [{i/'^)^vf]/e 



(7) 



to account for the friction of gas drag Vf acting on dust 
particles. As in Eq.®, the integral $ is a function of 
k. Our derivation of Eg. ([7]) was provided in the supple- 
mentary material of [20|. In this paper, we will carry out 
simulation tests to validate the use of Eq. ([T]) for a wide 
range of k. This validation test will be performed for the 
modest level of gas friction Vf vcl our experiment [20| . 

The TCAF measures the memory of transverse cur- 
rent, which refiects the decay of microscopic velocity 
shear. The shear decay can be caused by several mecha- 
nism in 2D dusty plasma clouds, such as Coulomb colli- 
sions amongst dust particles and the friction due to gas 
drag. We will study how gas friction affects the TCAF 
later. 



HI. SIMULATION METHODS 

In order to test the effects of gas friction, we will com- 
pare the results of two simulations: a Langevin MD sim- 
ulation with friction, and a frictionless equilibrium MD 
simulation. 

Our two simulation methods are the same in many re- 
spects. Both use a binary interparticle interaction with a 
Yukawa pair potential. In both simulations, particles are 
only allowed to move in a single 2D plane. Conditions re- 
mained steady during each simulation run. For both sim- 
ulations, the parameters we used were N = 4096 particles 
in a rectangular box with periodic boundary conditions. 
The box had sides 64.16 x 55.56. The integration time 
step was 0.019 w^^rf , and simulation data were recorded 
for a time duration of 68 000 after a steady state 
was reached. Both of our simulations were performed at 
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r = 68 and k = 0.5, which are the same values as in our 
experiment (20| . 

Our Langevin MD simulation takes into account the 
dissipation due to gas friction. The equation of motion 
that is integrated in the Langevin simulation is [36, 43, 



4*5 



mri = -V ^ 4>i.j ~ Vfinvi + C,i{t), 



(8) 



where vfmTi is a frictional drag and Ci(^) is a random 
force. There is no thermostat to adjust the tempera- 
ture; instead the temperature is established by choosing 
the magnitude of Q (t). Here, we chose the experimental 
value = 0.08 ujpd [20j. Note that this gas friction level 
is modest, i.e., the dust particle motion is underdamped, 
since i^/ <C ojpd- 

Our frictionless equilibrium MD simulation 43 , 5^, 13] 
has no gas friction in the equation of motion 



mr, 



(9) 



A Nose-Hoover thermostat is applied to maintain a de- 
sired temperature (sc " 



54| 



Trajectories ri{t) are found by integrating Eq. ([8|) or 
^ for all particles. An example is shown in Fig. 1 from 
the frictionless MD simulation. 



IV. RESULTS 
A. Hydrodynamic and viscoelastic regimes 

Comparing the results from the two simulations. Fig. 2, 
we can see how friction speeds the loss of memory of 
the system's microscopic shearing motion. The memory 
of the shearing motion is indicated by the decay of the 
TCAF. 

As expected (l^l , in the typical hydrodynamic limit of 
long length scales, as shown in Fig. 2(a), the TCAF is 
just a monotonic decay from unity to zero without any 
oscillations |23| . We find that at the same hydrodynamic 
length scale, the TCAF decays much faster with friction 
than without, indicating that in experimental dusty plas- 
mas gas friction plays an important role in shear decay 
in large length scales. 

When the wave number k is slightly larger, in the inter- 
mediate regime between the hydrodynamic and viscoelas- 
tic regimes, Fig. 2(b), the difference in TCAF between 
frictional and frictionless is smaller. The integral of the 
frictional TCAF is about a half of that for the friction- 
less TCAF, as seen in the inset of Fig. 2(b). This integral 
corresponds to as in Eq. © or Eq. ([7]). 

When the wave number k is even larger, in the vis- 
coelastic regime, the TCAF oscillates around zero af- 
ter its decay due to the elastic effects. Fig. 2(c). In 
this viscoelastic regime, there is little difference between 



the TCAF from the two simulations, indicating that at 
smaller length scales, gas friction does not contribute 
much to shear decay. The friction plays a larger role in 
TCAF at larger length scales than smaller length scales. 

The calculation of ry(fc) using Eq. ([6]) or ([7|) requires 
choosing an upper limit in the time integral of TCAF 
CT(fc,t). An infinite time is of course impractical for 
both experiments and simulations, so for a finite value 
we chose ti, the time of the first upward zero crossing 
of the TCAF as shown in Fig. 2(c). This choice is 
suitable for two reasons: first, it is sufficiently long to 
retain both viscous and elastic effects; second, we found 
that contributions to the integral after tj are negligible, 
for a TCAF that is not noisy. The calculation result 
for r]{k) is not very sensitive to the chosen upper limit. 
Extending the limit to a higher value would only cause a 
limited effect on the value of the integral. 



B. Validating the generalized ri{k) expression 



Results for the wave-number-dependent viscosity r](k) 
are presented in Fig. 3(b) and (c) for both simulations. 

We find an agreement in the values of r](k) for the 
frictionless and Langevin simulations. This agreement 
can be seen by comparing the circles in Fig. 3(b) for the 
frictionless simulation with Eq. ([6]), and the triangles in 
Fig. 3(c) for the Langevin simulation with Eq. 0. There 
is not only a qualitative agreement in the downward trend 
as the wave number k increases, but also a quantitative 
agreement. This quantitative agreement is most easily 
seen by fitting the calculated ri{k) to the Fade approx- 
imant of ^20„ , 39] and comparing the fit parameters, as 
indicated in Fig. 3 for the smooth curves. 

This agreement leads us to our first chief result: a 
validation of Eq. ^ for computing ri{k) in the presence of 
gas friction. Since the two simulations were performed for 
the same values of F and k, an agreement indicates that 
Eq. ([7]) is valid. If there had been a discrepancy between 
the circles in Fig. 3(b) and the triangles in Fig. 3(c), 
we would question whether Eq. ([7|) is valid. We gain 
confidence in the validity of Eq. ([7]) by the lack of any 
significant discrepancy in the two results. 

The importance of correcting for friction, in Eq. ([7]), is 
demonstrated in Fig. 3(c). If we use Eq. ([6]) instead, the 
presence of friction leads to an exaggerated value for i]{k), 
as seen by comparing the two sets of data in Fig. 3(c). 
This exaggeration is most extreme at small wave numbers 
(where the effect of friction is greatest, as we found in 
Sec. IV A for the TCAF). 
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C. Testing the Green-Kubo relation for static 
viscosity in the presence of friction 

To determine whether the Green-Kubo relation, 
Eq. (|4]), still provides an accurate calculation of static 
viscosity 77 of a 2D Yukawa liquid, in the presence of a 
modest level of gas friction, we performed a test of Eq. Q 
comparing rj computed from our frictional Langevin sim- 
ulation and our frictionless simulation. These results 
for the normalized kinematic static viscosity are Tj/p — 
(0.26 ± 0.02) a^uipd for the frictionless simulation, and 
rj/p — (0.27 ± 0.02) a^Wpd for the Langevin simulation 
with friction. These values are also shown in Fig. 3(b) 
and (c) as star symbols. Noting that these results are in 
agreement within the uncertainty, we conclude that the 
Green-Kubo relation remains accurate, at least with a 
modest level of gas friction, for a 2D Yukawa liquid at 
the value F = 68 and k — 0.5. 

A further confirmation of the accuracy of the Green- 
Kubo relation when used with modest levels of friction 
can be found by examining our ri{k) in Fig. 3(c). We note 
an agreement oi ri{k) as k ^ with ry from the Green- 
Kubo relation. This agreement is significant because ri{k) 
is computed from the TCAF, which is unrelated to the 
Green-Kubo relation used to compute rj. 

We can provide two intuitive suggestions to explain 
the accuracy of the Green-Kubo relation in the presence 
of a modest level of gas friction. First, we note that 
the gas friction that we have considered is so small that 
Vf/uipd < 0.1. This inequality demonstrates that fric- 
tional effects will in general be much smaller than effects 
arising from particle charge as measured by ujpd- Sec- 
ond, the TCAF in Fig. 2 showed us that gas friction has 
the least effect on motion at small length scales, and dy- 
namical information at these small length scales are also 
refiected in the Green-Kubo relation because it is based 
only on fiuctuations of individual particle motion. 

We cannot rule out the possibility that friction will af- 
fect the static viscosity computed using the Green-Kubo 
relation in other parameter regimes. In fact, for a 3D 
Yukawa Langevin simulation at a much lower F = 2, 
Ramazanov and Dzhumagulova found that 77 computed 
using the Green-Kubo relation diminishes as the friction 
was raised to a very high level (Hsj . 

IV. SUMMARY 

Motivated by experiments with 2D clouds of charged 
dust particles suspended in a plasma, we carried out two 
types of simulations, with and without gas friction. We 
validated the newly-introduced Eq. (O for calculating 
r]{k) as a measure of viscoelasticity, in the presence of 
gas friction. We also verified that the Green-Kubo rela- 
tion can accurately measure the static viscosity rj of the 
2D collection of charged dust particles even when they 



experience gas friction. The level of gas friction we con- 
sidered was at a low level Vf/ujpd < 0.1, and the coupling 
was moderate with F = 68 and k = 0.5, both as in our 
recent experiment [20| . 

This work was supported by NSF and NASA. 
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FIG. 1: (color online). Trajectories of particles from the fric- 
tionless MD simulation. Similar trajectories for the experi- 
ment and a Langevin MD simulation were reported in [20l |. 
All simulation results here are for F — 68, k — 0.5. 



FIG. 2: (color online). Transverse current autocorrelation 
functions (TCAF) from simulations in different regimes: (a) 
hydrodynamic, (b) intermediate between hydrodynamic and 
viscoelastic, and (c) viscoelastic. Friction plays a larger role 
at smaller k, i.e., longer length scales. After a long time, the 
TCAF always decays to zero. (For the small k case without 
friction (a), the TCAF approaches zero after a great time, 
longer than shown here.) The inset in (b) is the cumulative 
time integral of TCAF. 



FIG. 3: (color online). Wave- number-dependent viscosity 
ry(fc) from (a) the experiment, (b) the Langevin MD simula- 
tion, and (c) the frictionless MD simulation. The agreement 
of the smooth curves in (b) and (c) validates Eq. (O for calcu- 
lating ri{k) in the presence of gas friction. Using Eq. ([6]) in the 
presence of gas friction (circle data points in (c)) would lead 
to an exaggerated value for ri{k), especially at longer length 
scales. Also shown in (b) and (c) is the calculated static vis- 
cosity Tj based on the Green-Kubo relation as indicated with 
star symbols. Panel (a) is reprinted from [20l ]. 



